/**
 * Implement the half edges mesh structure.  
 * @author Songgang Xu (sxu at tamu dot edu)
 * @file mesh.h
 */

#ifndef XMESH_MANIFOLD_H
#define XMESH_MANIFOLD_H

#include <map>
#include <string>
#include <sstream>
#include <fstream>
#if defined(_MSC_VER)
#include <hash_map>
#elif defined(__GNUC__)
#include <ext/hash_map>
#endif

#include "primitive.h"

/**
 * Provide mesh interfaces for differential geometry. 
 * @file mesh.h
 */


template <typename Value> class HalfEdge ;
template <typename Value> class Edge;
template <typename Value> class Face;
template <typename Value> class Vertex; 
template <typename MeshPrimitive> class Mesh;

typedef std::pair<unsigned, unsigned> HashKey;
typedef HashKey HalfEdgeHashKey ;
typedef HashKey EdgeHashKey; 
class FaceHashKey
{
public:
	FaceHashKey(unsigned a, unsigned b, unsigned c)
	{
		idx1 = a; idx2 = b; idx3 = c;
	}
public:
	unsigned idx1, idx2, idx3 ;
};

#if defined(_MSC_VER)

template<>
struct stdext::hash_compare<HashKey > {
    // the mean number of elements per bucket
    static const size_t bucket_size = 4;
    // the minimum number of buckets
    static const size_t min_buckets = 65535*16;
    // hash function
    size_t operator()(const HashKey & key) const {
        return (key.first * key.second);
    }
    // compare function
    bool operator()(const HashKey & lhs, const HashKey & rhs) const {
        return lhs < rhs;
    }
};

/*
template<>
struct stdext::hash_compare<FaceHashKey > {
    // the mean number of elements per bucket
    static const size_t bucket_size = 4;
    // the minimum number of buckets
    static const size_t min_buckets = 65535*2;
    // hash function
    size_t operator()(const FaceHashKey & key) const {
		return key.idx1 * key.idx2 + key.idx3 ;
    }
    // compare function
    bool operator()(const FaceHashKey & lhs, const FaceHashKey & rhs) const {
		return (lhs.idx1 < rhs.idx1 || 
				(lhs.idx1 == rhs.idx1 && lhs.idx2 < rhs.idx2) ||
				(lhs.idx1 == rhs.idx1 && lhs.idx2 == rhs.idx2 && lhs.idx3 < rhs.idx3)) ;
    }
};
*/

#elif defined(__GNUC__)
struct EdgeHashOperator {
	size_t operator()(const EdgeHashKey & key) const {
        return (key.first * key.second);
    }
};
struct EdgeHashEqual {
	bool operator()(const EdgeHashKey & lhs, const EdgeHashKey & rhs) const {
		return lhs == rhs;
    }
};
/*
struct FaceHashOperator {
	size_t operator()(const FaceHashKey & key) const {
        return (key.idx1 * key.idx2 + key.idx3);
    }
};
struct FaceHashEqual {
	bool operator()(const FaceHashKey & lhs, const FaceHashKey & rhs) const {
		return ( lhs.idx1 == lhs.idx1 
				 && lhs.idx2 == lhs.idx2 
				 && lhs.idx3 == lhs.idx3 );
    }
};
*/
#endif
template <typename Value>
inline HalfEdgeHashKey halfEdgeHashKey (HalfEdge<Value> * he)
{
	unsigned id1 = he->from_v->idx ;
	unsigned id2 = he->next->from_v->idx; 
	return std::make_pair (id1, id2);
}
inline EdgeHashKey edgeHashKey (unsigned id1, unsigned id2)
{
	return (id1 < id2 ? std::make_pair(id1, id2) : std::make_pair(id2, id1)) ;
}
template <typename Value>
inline EdgeHashKey edgeHashKey(Edge<Value> * edge)
{
	if (edge == NULL) 
	{
		std::cout<<"INVALID EDGE"<<std::endl;
	}
	unsigned id1 = edge->he->from_v->idx;
	unsigned id2 = edge->he->next->from_v->idx ;
	return edgeHashKey(id1, id2);
}

template <typename Value>
struct MeshPrimitive
{
	typedef Value            ValueType ; 
	typedef Face <Value>     FaceType ;
	typedef Vertex <Value>   VertexType  ;
	typedef HalfEdge <Value> HalfEdgeType ; 
	typedef Edge <Value>     EdgeType ;
	typedef Vector<Value>    VectorType ;
	typedef VectorType       PointType ;

};

/**
 * HalfEdge Class. 
 */
template <typename Value>
class HalfEdge 
{
protected:
	typedef typename MeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename MeshPrimitive<Value>::FaceType     FaceType;
	typedef typename MeshPrimitive<Value>::VertexType   VertexType;
	typedef typename MeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename MeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename MeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename MeshPrimitive<Value>::PointType    PointType ;

public:
	/**
	 * HalfEdge Constructor.
	 */
	HalfEdge (): next(NULL), prev(NULL), flip(NULL), face(NULL), from_v(NULL), edge(NULL){}
	HalfEdgeType * next, * prev, * flip;
	FaceType * face;
	VertexType * from_v ;
	EdgeType * edge; 

};



/**
 * Vertex Class. 
 */
template<typename Value>
class Vertex
{
protected:
	typedef typename MeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename MeshPrimitive<Value>::FaceType     FaceType;
	typedef typename MeshPrimitive<Value>::VertexType   VertexType;
	typedef typename MeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename MeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename MeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename MeshPrimitive<Value>::PointType    PointType ;
	//typedef Vector<Value> PointType ;

public:

	/**
	 * Constructor. 
	 * @param _coord Coordinates for Vertex.
	 * @param _texcoord Texture Coordinates for Vertex .
	 * @param _idx Index for Vertex. 
	 */
    Vertex(const PointType & _coord = PointType (ValueType(0), ValueType(0), ValueType(0)),  unsigned idx = -1, const PointType & _texcoord = PointType (ValueType(0), ValueType(0))) 
		: m_idx(idx), m_coord(_coord), m_texcoord(_texcoord) 
	{
		//this->m_local_area = ValueType(0); 
	}
	
public:
	
	/**
	 * Calculate the area of local averaging region. 
	 */
	const ValueType get_local_area ()  const 
	{
		ValueType area = ValueType(0);
		PointType center (ValueType(0), ValueType(0), ValueType(0));
		unsigned count = 0; 
		for (VVIter iter = vv_begin(); iter != vv_end(); ++iter)
		{
			//Fix me later. I will define a new class Polygon. 
			center += iter->m_coord ;
			++count ;
		}
		center /= ValueType(count); 


		
		return this->m_local_area ; 
	}
	/**
	 * Vertex Face Iterator. 
	 */
	class VFIter 
	{
	public:
		VFIter (VertexType * v)
		{
			he = v->he; 
		};
		VFIter (HalfEdgeType * edge)
		{
			he = edge; 
			begin = edge; 
		}
		VFIter & operator ++()		
		{
			if (this->he->flip == NULL)	
				this->he = begin ;
			else
				this->he = he->flip->next ;
			return (*this) ;
		}
		FaceType * operator * ()
		{
			return (this->he->face);
		}
		bool operator ==(const VFIter & other)
		{
			return(he == other.he);
		}

		bool operator !=(const VFIter & other)
		{
			return(he != other.he);
		}

	private :
		HalfEdgeType * he; 
		HalfEdgeType * begin ;
	};
	//This two iterator will deal with the mesh has boundary
	VFIter vfbegin_safe()
	{
		//Check wheter the vertex is on boundary or not
		HalfEdgeType * iter  = this->he;
		HalfEdgeType * iterEnd  = this->he;
		do
		{
			if (iter->prev->flip == NULL)
			{
				bound = iter ;
				return VFIter (bound) ;
			}
			iter = iter->prev->flip ;
		}while (iter != iterEnd) ;
		bound = iterEnd ;
		return VFIter (bound) ;
	}
	VFIter vfend_safe()
	{
		return VFIter (bound);
	}
	VFIter vfbegin ()
	{
		return VFIter (this) ;
	}
	VFIter vfend ()
	{
		return VFIter (this->he->prev->flip);
	}
public:
	/**
	 * Vertex Vertex Iterator. 
	 */
	class VVIter 
	{
	public:
		VVIter (VertexType * v)
		{
			he = v->he; 
			begin = he; 
		};
		VVIter (HalfEdgeType * edge)
		{
			he = edge; 
			begin = he ;
		}
		VVIter & operator ++()
		{
			this->he = he->flip->next ;
			if (he == begin) he = NULL ;
			return (*this) ;
		}
		VertexType * operator*()
		{
			return (this->he->flip->from_v);
		}
		bool operator ==(const VVIter & other)
		{
			return(he == other.he);
		}

		bool operator!=(const VVIter & other)
		{
			return(he != other.he);
		}
		const  HalfEdgeType * get_he () const {return this->he ;}
	private :
		HalfEdgeType * he;
		HalfEdgeType * begin ;
	};
public:
	VVIter vv_begin ()
	{
		return VVIter (this) ;
	}
	VVIter vv_end ()
	{
		return VVIter ((HalfEdgeType * ) (NULL));
	}

public: 
	int m_idx ;                // Vertex Index 
	PointType m_coord ;        // Coordinates for vertex
	PointType m_texcoord;      // Texture coordinates for vertex
	PointType m_normal ;       // Vertex normal    
	HalfEdgeType * he;         // HalfEdge structure. 
	HalfEdgeType * bound ;     //Used for mesh with boundary

//	static void * mp_mesh ; 
};



/**
 * Edge class.
 *
 */
template <typename Value>
class Edge 
{
protected:
	typedef typename MeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename MeshPrimitive<Value>::FaceType     FaceType;
	typedef typename MeshPrimitive<Value>::VertexType   VertexType;
	typedef typename MeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename MeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename MeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename MeshPrimitive<Value>::PointType    PointType ;

public :
	HalfEdgeType * he ; 
	friend std::ostream & operator << (std::ostream & ostr, EdgeType * e)
	{
		ostr <<'('<<e->he->from_v->idx<<','; 
		ostr <<e->he->flip->from_v->idx<<')';
		return ostr;
	}
	Value cost; 
//	static void * mp_mesh; 
};



/**
 * Face class.
 *
 */
template <typename Value>
class Face 
{
protected:
	typedef typename MeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename MeshPrimitive<Value>::FaceType     FaceType;
	typedef typename MeshPrimitive<Value>::VertexType   VertexType;
	typedef typename MeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename MeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename MeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename MeshPrimitive<Value>::PointType    PointType ;


		
public :
	/**
	 * Face HalfEdge Iterator
	 */
	class FHeIter
	{
	public:
		FHeIter (FaceType * face)
		{
			he = face->he; 
			begin = he; 
		}
		FHeIter (HalfEdgeType * edge)
		{
			he = edge; 
			begin  = he; 
		}
		FHeIter & operator ++()
		{
			this->he = he->next ;
			if (he == begin) he = NULL ;
			return (*this) ;
		}
		HalfEdgeType * operator*()
		{
			return (this->he);
		}
		bool operator ==(const FHeIter & other)
		{
			return(he == other.he);
		}
		bool operator != (const FHeIter & other)
		{
			return(he != other.he);
		}

	private:
		HalfEdgeType * he    ; 
		HalfEdgeType * begin ;
	}; 

	FHeIter fhe_begin ()
	{
		return FHeIter(this);
	}
	FHeIter fhe_end ()
	{
		return FHeIter((HalfEdgeType * ) (NULL));
	}
public :
	/**
	 * Face Vertex Iterator
	 */
	class FVIter
	{
	public:
		FVIter (const FaceType * face)
		{
			he = face->he; 
			begin = he; 
		}
		FVIter (HalfEdgeType * edge)
		{
			he = edge; 
			begin  = he; 
		}
		FVIter & operator ++()
		{
			this->he = he->next ;
			if (he == begin) he = NULL ;
			return (*this) ;
		}
		VertexType * operator*()
		{
			return (this->he->from_v);
		}
		bool operator ==(const FVIter & other)
		{
			return(he == other.he);
		}
		bool operator != (const FVIter & other)
		{
			return(he != other.he);
		}
		HalfEdgeType * get_he () const
		{
			return this->he; 
		}
	private:
		HalfEdgeType * he    ; 
		HalfEdgeType * begin ;
	}; 

	FVIter fv_begin () const 
	{
		return FVIter(this);
	}
	FVIter fv_end () const 
	{
		return FVIter((HalfEdgeType * ) (NULL));
	}
public:
	Face () : m_idx (-1)
	{

	}

	/**
	 * Get the normal of the face. 
	 * @return normal of the face. 
	 */
	const VectorType & normal () const  
	{
		return this->normalVector ; 
	}

    /**
	 * Get Vertex coordinates on the face in counter-clockwise order. 
	 * @param v_set Vertex coordinates index set. 
	 */
	void Face2VertexCoord (std::vector<PointType *> & v_set)
	{
		v_set.clear ();
		for (FVIter iter = fv_begin(); iter != fv_end(); ++iter)
		{
			v_set.push_back( &((*iter) ->m_coord));
		}
	}

	/**
	 * Get Vertex Pointer the face in counter-clockwise order. 
	 * @param v_set Vertex index set. 
	 */
	inline void Face2Vertex (std::vector<VertexType *> & v_set)
	{
		v_set.clear ();
		for (FVIter iter = this->fv_begin(); iter != this->fv_end(); ++iter)
		{
			v_set.push_back((VertexType *) (*iter));
		}
	}


	/**
	 * Get Vertex Index on the face in counter-clockwise order. 
	 * @param v_set Vertex index set. 
	 */
	inline void Face2VertexIdx (std::vector<int> & v_set)
	{
		v_set.clear ();
		for (FVIter iter = fv_begin(); iter != fv_end(); ++iter)
		{
			v_set.push_back( (*iter)->m_idx);
		}
	}

    /**
	 * Get Vertex texture coordinates on the face in counter-clockwise order. 
	 * @param v_set Vertex texturecoordinates index set. 
	 */
	inline void Face2TexCoord (std::vector<PointType *> & v_set)
	{
		v_set.clear ();
		for (FVIter iter = fv_begin(); iter != fv_end(); ++iter)
		{
			v_set.push_back( &((*iter) ->m_texcoord));
		}
	}

	inline void set_normal (const VectorType & _n)
	{
		this->normalVector = _n ;
	}
	/**
	 * Calculate the face nroaml ; 
	 */
	inline void compute_normal ()
	{
		std::vector<PointType * > p ;
		Face2VertexCoord(p);
		VectorType p1 = (*p[1]) -  (*p[0]);
		VectorType p2 = (*p[2]) -  (*p[0]);
		normalVector = (p1^p2).normalize();			
	}

	/**
	 * This function could be inherited to construct more data sturcture on faces. 
	 */
	inline virtual void construct()
	{
		compute_normal();
	}

	/**
	 * This function is called to update some data structructures after some operations. 
	 * For example the transplate operation. 
	 */
	inline virtual void update ()
	{
		//this->construct();
	}

	virtual ~Face() {}

public: 

	int m_idx ;
	HalfEdgeType * he;
	VectorType normalVector ;
//	static void * mp_mesh ; 

} ;

/**
 * A shared point points to the mesh. 
 */
/*
template <typename Value> void * Edge<Value>::mp_mesh = NULL ;
template <typename Value> void * Vertex<Value>::mp_mesh = NULL ;
template <typename Value> void * HalfEdge<Value>::mp_mesh = NULL ;
template <typename Value> void * Face<Value>::mp_mesh = NULL ;
*/
/**
 * Mesh class. 
 */
template <typename MeshPrimitive>
class Mesh
{
protected:
	typedef typename MeshPrimitive::ValueType    ValueType ; 
	typedef typename MeshPrimitive::FaceType     FaceType ;
	typedef typename MeshPrimitive::VertexType   VertexType ;
	typedef typename MeshPrimitive::HalfEdgeType HalfEdgeType ; 
	typedef typename MeshPrimitive::EdgeType     EdgeType ; 
	typedef typename MeshPrimitive::VectorType   VectorType ;
	typedef typename MeshPrimitive::PointType    PointType ;


#if defined(_MSC_VER)
	typedef stdext::hash_map <unsigned, VertexType * >  VertexContainer ;
	typedef stdext::hash_map <EdgeHashKey, EdgeType * > EdgeContainer ;
	typedef stdext::hash_map <HashKey, HalfEdgeType * > HalfEdgeContainer ;
	typedef std::vector      <FaceType * >              FaceContainer ;
	
#elif defined(__GNUC__)
	typedef __gnu_cxx::hash_map <unsigned, VertexType * >                                   VertexContainer ;
	typedef __gnu_cxx::hash_map <EdgeHashKey, EdgeType * , EdgeHashOperator, EdgeHashEqual> EdgeContainer ;
	typedef __gnu_cxx::hash_map <HashKey, HalfEdgeType * , EdgeHashOperator, EdgeHashEqual> HalfEdgeContainer ;
	typedef std::vector    <FaceType *>                                                     FaceContainer ;
#endif

public :
	typedef typename VertexContainer::iterator VertexIterator ;
	typedef typename EdgeContainer::iterator EdgeIterator ;
	typedef typename FaceContainer::iterator FaceIterator ;
	typedef typename HalfEdgeContainer::iterator HalfEdgeIterator ;
	typedef typename FaceType::FVIter            FVIter ;

public:
	Mesh() {}
	virtual ~Mesh() {this->finalize ();}
	friend std::istream & operator >> (std::istream & istr, Mesh & mesh)
	{
		return mesh.OBJ_loader(istr );
		//return mesh.OFF_loader(istr ); 
	}
	/**
	 * OBJ File reader .
	 */
	virtual void OBJ_reader (std::istream & istr)
	{
		this->OBJ_loader (istr );
	}

	/**
	 * OBJ File Writer. 
	 */
	virtual void OBJ_writer (std::ostream & ostr)
	{
		for (unsigned i = 0 ; i < this->vertices.size(); ++i)
		{
			PointType p = vertices[i]->m_coord ;
			ostr<<"v "<<p[0]<<' '<<p[1]<<' '<<p[2] <<std::endl; 
		}
		for (unsigned i = 0; i < this->faces.size(); ++i)
		{
			ostr<<"f"; 
			for (FVIter iter = this->faces[i]->fv_begin(); iter != this->faces[i]->fv_end(); ++iter)
			{
				ostr<<' '<<((*iter)->m_idx+1) ;
			}
			ostr<<std::endl;
		}
	}

	/**
	 * OFF File reader. 
	 */
	virtual void OFF_reader (std::istream & istr)
	{
		this->OFF_loader (istr);
	}

	/**
	 * OFF Writer ()
	 */
	virtual void OFF_writer ()
	{
		throw ("Incomplete function. ");
	}


public:
	/**
	 * Rotate the model 
	 * Use rotation matrix. 
	 */
	void rotateX ( const ValueType & rot)
	{
		for (VertexIterator v_iter = vertices.begin(); v_iter != vertices.end(); ++ v_iter)
		{
			v_iter->second->m_coord = VectorRotateX3D (v_iter->second->m_coord, rot ) ; 
		}
	}
	/**
	 * Translate the model ;
	 */
	void translate (const VectorType & _v)
	{
		for (VertexIterator v_iter = vertices.begin(); v_iter != vertices.end(); ++ v_iter)
		{
			v_iter->second->m_coord += _v; 
		}
		for (FaceIterator f_iter = faces.begin(); f_iter != faces.end(); ++ f_iter )
		{
			(*f_iter)->update (); 
		}
	}
	
public:
	/**
	 * Interface for calculating curvature at position _v. 
	 * @param _v Vector indicates position coordinates. 
	 */
	virtual ValueType curvature (const VectorType & _v) 
	{
		throw ("Incomplete function");
		return ValueType(0);
	} 

	/**
	 * Interface for calculating gradient at position _v. 
	 * @param _v Vector indicates position coordinates. 
	 */
	virtual VectorType gradient (const VectorType & _v) 
	{
		throw ("Incomplete function");
		return VectorType (ValueType(0), ValueType(0), ValueType(0));
	}

	/**
	 * Interface for calculating tensor at position _v. 
	 * @param _v Vector indicates position coordinates. 
	 */
	virtual VectorType tensor (const VectorType & _v) 
	{
		throw ("Incomplete function");
		return VectorType (ValueType(0), ValueType(0), ValueType(0));
	}

	/**
	 * Interface for calculating local area at position _v. 
	 * @param _v Vector indicates position coordinates. 
	 */
	virtual VectorType area (const VectorType & _v) 
	{
		throw ("Incomplete function");
		return VectorType (ValueType(0), ValueType(0), ValueType(0));
	}

private:
	istream & OFF_loader (std::istream & istr)
	{
		std::string line ;
		std::getline(istr, line);
		if (line != "OFF") return istr ;

		unsigned vertex_count = 0;
		unsigned v_num , face_num, edge_num; 
		VertexType * v; 
		while (std::getline(istr, line))
		{
			if (line.size() >0 && line[0] == '#') continue; 
			std::stringstream sstr(line.c_str());
			sstr >> v_num >> face_num >> edge_num;
			
			for (unsigned i = 0; i < v_num ; ++i)
			{
				v = new VertexType ;
				istr>>v->m_coord ;
				v->he = NULL ;
				this->vertices[vertex_count] = v;
				++ vertex_count ;
			}
			for (unsigned i = 0; i < face_num ; ++i)
			{
				unsigned numOfVertex; 
				std::vector <unsigned > v_id ;
				istr>>numOfVertex ;
				for (unsigned j = 0; j < numOfVertex ; ++j)
				{
					unsigned v_idx; 
					istr>> v_idx; 
					v_id.push_back(v_idx);
				}
				std::getline(istr, line);
				construct_face(v_id); 
			}
			break; 
		}

		/** Initialization of the shared pointer*/
		construct();
		return istr; 
	}
	istream & OBJ_loader (std::istream & istr)
	{
		std::string line;
		std::string tag; 
		std::string slash;
		VertexType * v; 
		unsigned count = 0; 
		unsigned tex_count = 0;

		while (std::getline(istr, line))
		{
			std::stringstream sstr(line.c_str());
			sstr >> tag;
			if (tag == "v")
			{	
				v = new VertexType;
				sstr >> v->m_coord ;
				v->he = NULL; 
				v->m_idx = count ;
				this->vertices[count]=v;
				++ count ;
			}
			else if (tag == "vt")
			{
				VectorType tex_coord(2);
				sstr>>tex_coord ;
				this->vertices[tex_count]->m_texcoord = tex_coord ;
				++ tex_count ;
			}
			else if (tag == "f")
			{
				std::vector <unsigned > v_id ;

                while(sstr) {
                    std::string s ;
                    sstr >> s ;
                    if(s.length() > 0) {
                        std::stringstream v_input(s.c_str()) ;
                        int index ;
                        v_input >> index ;
						v_id.push_back(index - 1) ;
                        char c ;
                        v_input >> c ;
                        if(c == '/') {
                            v_input >> index ;
                        }
                    }
                }
				this->construct_face(v_id);
			}
		}
		construct();
		return istr ;
	}
	/**
	 * Post constructraction after the mesh is constructed. 
	 */
	virtual void construct () {}

protected: 

	/**
	 * Construct the mesh face.
	 * This function could be inherited to construct more data structures. 
	 */
	FaceType * construct_face (const std::vector <unsigned > & v_id) 
	{
		unsigned v_size = v_id.size(); 
	 	FaceType * p_face = new FaceType; 
		std::vector<HalfEdgeType * > v_he; 
		for (unsigned i = 0; i < v_size; ++ i)
		{
			v_he.push_back(new HalfEdgeType);
		}

		for (unsigned i = 0; i < v_size; ++i)
		{
			v_he[i]->face = p_face ;
			v_he[i]->prev = v_he[(i + v_size - 1) % v_size] ;
			v_he[i]->next = v_he[(i + 1) % v_size] ;
			
			v_he[i]->from_v = vertices[v_id [i]];
			vertices[v_id[i]]->he = v_he[i] ;

			hash_he [HashKey (v_id[i], v_id[(i + 1 ) % v_size])] = v_he[i] ;
		}

		p_face->he = v_he[0] ;
		p_face->m_idx =  this->faces.size();
		this->faces.push_back(p_face);


		// Update the flips and Edges 
		for (unsigned i = 0; i < v_size; ++i)
		{
			unsigned id1 = v_id[i];
			unsigned id2 = v_id[ (i + 1) % v_size] ;

			HalfEdgeIterator iter ;

			iter = hash_he.find(HashKey(id2, id1)) ;
			if (hash_he.end() == iter)
			{
				v_he[i]->flip = NULL;
				EdgeType * edge = new EdgeType;
				edge->he = v_he[i] ;
				v_he[i]->edge = edge;
				if (id1 < id2) edges[HashKey(id1, id2)] = edge;
				else edges[HashKey(id2, id1)] = edge;
			}
			else
			{
				v_he[i]->flip = (*iter).second;
				(*iter).second->flip = v_he[i] ;
				v_he[i]->edge = (*iter).second->edge ;
			}
			
		}

		p_face->construct(); 
		return p_face; 
	}

private :
	void finalize ()
	{
		for (VertexIterator it = vertices.begin();
			 it != vertices.end(); ++it)
		{
			if (it->second) delete (it->second);
			it->second = NULL; 			
		}
		std::cerr<<"Remove vertices"<<std::endl;

		for (FaceIterator it = faces.begin();
			 it != faces.end(); ++it)
		{
			if (*it) delete (*it);
		    (*it) = NULL; 			
		}
		std::cerr<<"Remove faces"<<std::endl;
		
		for (EdgeIterator it = edges.begin(); 
			 it !=edges.end(); ++it)
		{
			if (it->second) delete it->second; 
			it->second = NULL; 
		}
		std::cerr<<"Remove edges"<<std::endl;
		for (HalfEdgeIterator it = hash_he.begin(); 
			 it != hash_he.end(); ++it)
		{
			if (it->second) delete it->second; 
			it->second = NULL; 
		}
		std::cerr<<"Remove Halfedges"<<std::endl;
	}

public:
	VertexContainer   vertices;
	EdgeContainer     edges; 
	FaceContainer     faces; 
	HalfEdgeContainer hash_he; 

};

/**
 * @page tutorial_mesh_construction Mesh Construction
 * @section tutorial_mesh_vv_iter Iterator for Vertex. 
 * In this example, we show how to use Vertex Vertex Iterator. The program outputs the neighbors of each vertex. 
 * @section tutorial_mesh_vv_iter_example Example
 * Example usage:
 * @code
#include <iostream>
#include <cmath>
#include "matrix.h"
#include "matrix_util.h"
#include "mesh.h"


typedef typename MeshPrimitive<float>::PointType     PointType  ;
typedef typename MeshPrimitive<float>::VectorType    VectorType ;
typedef typename MeshPrimitive<float>::FaceType      FaceType   ;
typedef typename MeshPrimitive<float>::VertexType    VertexType   ;
typedef Mesh <MeshPrimitive<float> >                 MeshType ;

int main(int argc, char ** argv)
{
        MeshType mesh ;
        try
        {

                std::ifstream ifs (argv[1]);
                ifs >> (mesh) ;
		
		for (MeshType::VertexIterator iter = mesh.vertices.begin(); iter != mesh.vertices.end(); ++iter)
		{
			VertexType * v = iter->second ;
			std::cout<<v->m_idx<<" : "; 
			for (VertexType::VVIter vv_iter = v->vv_begin(); vv_iter != v->vv_end(); ++vv_iter)
			{
				std::cout<<(*vv_iter)->m_idx<<' ';
			}
			std::cout<<std::endl;
		}
        }
        catch (const char * msg)
        {
                std::cerr<<msg<<std::endl;      
        }
}
 * @endcode
 */

/**
 * @page tutorial_mesh_construction Mesh Construction
 *
 * @section tutorial_mesh_con_intr Introduction
 * In this example, we show how to extend the Mesh Primitives by adding customized attributes and member functions. 
 @section tutorial_mesh_con_example Example
 * Example usage:
 @code
#include <iostream>
#include <cmath>
#include "matrix.h"
#include "matrix_util.h"
#include "mesh.h"

template <typename Value>
class SubHalfEdge : public HalfEdge<Value>
{
public :
	SubHalfEdge () : weight (0) , HalfEdge <Value>() {}
public :
	unsigned weight; 
	unsigned get_weight() {return this->weight ;}
};

template <typename Value>
struct SubMeshPrimitive : public MeshPrimitive <Value>
{
//	typedef Value            ValueType ; 
//	typedef Face <Value>     FaceType ;
//	typedef Vertex <Value>   VertexType  ;
	typedef SubHalfEdge <Value> HalfEdgeType ; 
//	typedef Edge <Value>     EdgeType ;
//	typedef Vector<Value>    VectorType ;
//	typedef VectorType       PointType ;

};


typedef typename SubMeshPrimitive<float>::PointType     PointType  ;
typedef typename SubMeshPrimitive<float>::VectorType    VectorType ;
typedef typename SubMeshPrimitive<float>::FaceType      FaceType   ;
typedef Mesh <SubMeshPrimitive<float> >                 MeshType ;

int main(int argc, char ** argv)
{
	MeshType mesh ;
	try
	{

		std::ifstream ifs (argv[1]);
		ifs >> (mesh) ;
	
		for (MeshType::HalfEdgeIterator iter = mesh.hash_he.begin(); iter != mesh.hash_he.end(); ++iter)
		{
			std::cout<<((*iter).second)->get_weight()<<' ';
		}
	}
	catch (const char * msg)
	{
		std::cerr<<msg<<std::endl;	
	}
}
 * @endcode
 * Text to immediately follow the code block.
 */

#endif
